import numpy as np

class Extract_2d(object):
    '''提取二维特征参数'''

    def __init__(self, z_data, x_coord):
        
        self.z_data = z_data		#三维数据
        self.x_coord = x_coord		#计算的行或列的序号
        self.num = z_data.shape[0]	#某行或列数据的个数
        
    #Rp表示..
    def extractRp(self,T):
        Rp=max(T)
        return Rp

    def extractRv(self,T):
        Rv1=min(T)
        Rv=abs(Rv1)
        return Rv

    def extractRz(self,T):
        Rp=max(T)
        Rv=min(T)
        Rz=Rp-Rv
        return Rz

    def extractRa(self,T):
        sum=0
        for i in range(T.shape[0]):
            sum = sum + abs(T[i])
            Ra=sum/T.shape[0]
        return Ra

    def extractRq(self,T):
        sum=0
        for i in range(T.shape[0]):
            sum = sum + T[i]*T[i]
            Rq=np.sqrt(sum/T.shape[0])
        return Rq
    def findmax(self,k,m,T):
        max=0
        for i in range(k,1,m):
            if max<abs(T[i]):
                max=abs(T[i])
        return max


    def extractRc(self,T):
        k=0
        m=0
        M=[]
        for i in range(len(T)-1):
            if T[i]*T[i+1]<=0:
                k=m
                m=i
                M.append(self.findmax(k,m,T))
        M.append(self.findmax(m,599,T))
        lenth=len(M)
        c=0
        if lenth%2==1:
            for i in range(lenth-1):
                c=c+M[i]
            Rc=c/(lenth-1)*2
        elif lenth%2==0:
            for i in range(lenth):
                c=c+M[i]
            Rc=c/(lenth/2)
        else:
            Rc=0
        return Rc
    


    #待解决零点的求解
    def extractRsm(self,T):
        mid = np.sum(T) / T.shape[0]
        p=float(mid+0.01)
        q=float(mid-0.01)
        m = np.zeros([T.shape[0]])
        k=0
        for j in range(T.shape[0]):
            if q<=float(T[j])<=p:
                m[k] = j#j为符合条件的点的位置
                k=k+1
        d = np.zeros([T.shape[0]])
        s=np.argmax(m)
        s=s+1
        l=m[s-1]-m[0]
        Rsm=2*l/s
        Rsm = Rsm/T.shape[0]#像素点个数转为长度
        return Rsm

    def extractRmrc(self,T,c=0.0015):
        sum=0
        for i in range(T.shape[0]):
            if T[i]>=c:
                sum=sum+1
        Rmr=sum/T.shape[0]
        return Rmr

    def extractRdq(self,T):
        sum=0
        for i in range(3,T.shape[0]-3):
            sum=sum+(1/60*(T[i+3]-9*T[i+2]+45*T[i+1]-45*T[i-1]+9*T[i-2]-T[i-3]))**2
        Rdq=np.sqrt(sum/594)
        return Rdq

    def extractrRsk(self,T, Rq):
        sum = 0
        for j in range(T.shape[0]):
            sum = sum + T[j]**3
        Rsk = sum/(T.shape[0]*Rq**3)
        return Rsk

    def extractrRku(self,T,Rq):
        sum = 0
        for j in range(T.shape[0]):
            sum = sum + T[j]**4
        Rku = sum/(T.shape[0]*Rq**4)
        return Rku


    def findhight(self,x,y,mr):
        for i in range(50):
            if y[i]<mr:
                hight=(x[i]+x[i-1])/2
                break
        return hight

    def extractRdc(self,T, mr1, mr2):  #要求mr1<mr2

        x = np.linspace(np.min(T), np.max(T), 50)

        y=[0 for i in range(50)]
        for i in range(50):
            y[i]=self.extractRmrc(T,x[i])

        for i in range(50):
            if y[i]<mr1:
                hight1=(x[i]+x[i-1])/2
                break
        for j in range(50):
            if y[j]<mr2:
                hight2=(x[j]+x[j-1])/2
                break
        Rdc=hight1 - hight2
        return Rdc

    def extractmidline(self,T):
        avg = np.sum(T) / T.shape[0]
        '''
        T = np.zeros(T.shape[0])
        avg = 0
        for i in range(T.shape[0]):
            T[i] = (float(T[i]))/T.shape[0]
            avg = avg + T[i]
        '''
        return avg


    def calc_2d_args(self):

        Rp = self.extractRp(self.z_data[self.x_coord])
        Rv = self.extractRv(self.z_data[self.x_coord])
        Rz = self.extractRz(self.z_data[self.x_coord])
        Ra = self.extractRa(self.z_data[self.x_coord])
        Rq = self.extractRq(self.z_data[self.x_coord])
        Rc = self.extractRc(self.z_data[self.x_coord])
        Rsm = self.extractRsm(self.z_data[self.x_coord])
        Rmrc = self.extractRmrc(self.z_data[self.x_coord])
        Rsk = self.extractrRsk(self.z_data[self.x_coord], Rq)
        Rku = self.extractrRku(self.z_data[self.x_coord], Rq)
        Rdq = self.extractRdq(self.z_data[self.x_coord])
        Rdc = self.extractRdc(self.z_data[self.x_coord], 0.1, 0.8)
        
        args_2d=(
                        '%.3f' % Rp,
                        '%.3f' % Rv,
                        '%.3f' % Rz,
                        '%.3f' % Rc,
                        '%.3f' % Ra,
                        '%.3f' % Rq,
                        '%.3f' % Rsm,
                        '%.3f' % Rmrc,
                        '%.3f' % Rsk,
                        '%.3f' % Rku,
                        '%.3f' % Rdq,
                        '%.3f' % Rdc)
        return args_2d




class Extract_3d(object):
    '''提取三维特征参数'''

    def __init__(self, z_data):
        self.z = z_data
        self.num = self.z.size
        
    def calc_sa(self):
        '''算数平均高度Sa'''
        return np.divide(np.sum(np.abs(self.z)), self.num)

    def calc_sz(self):
        '''最大高度'''
        return np.max(self.z) - np.min(self.z)

    def calc_sq(self):
        '''尺寸限定表面的均方根高度'''
        return np.sqrt(np.divide(np.sum(self.z ** 2), self.num))
        
    def calc_ssk(self):
        '''尺寸限定表面的偏斜度'''
        sq = self.calc_sq()
        return np.divide(np.sum(self.z ** 3), self.num* sq ** 3)
        

    def calc_sku(self):
        '''尺寸限定表面的陡峭度'''
        sq = self.calc_sq()
        return np.divide(np.sum(self.z ** 4), self.num* sq ** 4)

    def calc_sp(self):
        '''最大峰高Sp'''
        return np.max(self.z)


    def calc_sv(self):
        '''最大谷深Sv'''
        return np.min(self.z)

    def calc_smr(self, *c):
        '''
        尺寸限定表面的支撑面积率

        parameters
        ----------
        c : list, tuple, or 以','分隔的高度
        
        returns
        -------
        smr : list
        '''
        #return np.divide(self.z[self.z >= c].size, self.num)
        smr = [np.divide(self.z[self.z >= i].size, self.num) for i in c]
        return smr



    def calc_vvv_vvc_vmp_vmc(self, p1=0.1, p2=0.8):
        '''
        峰部的实体体积vmp: 表示负载面积率p1(10%)的实体体积。
        谷部的空隙容积vvv：表示负载面积率p2(80%)的谷部的空隙容积。
        '''
        z_min, z_max = np.min(self.z), np.max(self.z)
        
        def fun_smr(c, p):
            return np.divide(self.z[self.z >= c].size, self.num) - p
        
        def binary_search(fun, p, start=z_min, end=z_max, precision=0.001):
            '''
            二分法查找smr=p时对应的高度值c

            parameters
            ----------
            fun : 函数名    p : smr的值     precision : 精度

            returns
            -------
            start : smr=p对应的p值

            '''
            fa = fun(start, p)
            fb = fun(end, p)
            nums = 0
            if abs(fa) < 0.001:     # 判断区间边界值是否符合条件，近似零视为等于0
                return start
            if abs(fb) < 0.001:
                return end
            while abs(end - start) > precision and nums < 20:
                mid = (start + end) / 2.0
                fmid = fun(mid, p)
                if fmid == 0:
                    return mid
                if fa * fmid < 0:
                    end = mid
                if fmid * fb < 0:
                    start = mid
                nums += 1
            return start
             
        c_vmp = binary_search(fun_smr, p1)  # smr=0.1时对应的高度
        c_vvv = binary_search(fun_smr, p2)  # smr=0.8时对应的高度
        
        # 积分法求面积vmp,vvv,vvc,vmc,步长取0.01
        step = 0.01
        x_vmp = list(np.arange(c_vmp, z_max, step))  #step=0.01
        x_vvv = list(np.arange(z_min, c_vvv, step))
        x_vmc = list(np.arange(c_vvv, c_vmp, step))
        
        vmp = np.sum((np.array(self.calc_smr(*x_vmp))) * step)
        vvv = np.sum((1 - (np.array(self.calc_smr(*x_vvv)))) * step)
        vmc = np.sum((np.array(self.calc_smr(*x_vmc))) * step)
        vvc = (c_vmp - c_vvv) * 1 - vmc

        return vvv, vvc, vmp, vmc
        
    def calc_sdq(self):
        
        dx = (self.z[2:, 1:-1] - self.z[:-2, 1:-1]) / 2.0
        dy = (self.z[1:-1, 2:] - self.z[1:-1,:-2]) / 2.0
        return np.sqrt(np.sum(dx ** 2 + dy ** 2) / self.num)
        

    def calc_sdr(self):
        dx = (self.z[2:, 1:-1] - self.z[:-2, 1:-1]) / 2.0
        dy = (self.z[1:-1, 2:] - self.z[1:-1,:-2]) / 2.0
        return np.divide(np.sum(np.sqrt(1 + dx ** 2 + dy ** 2) - 1), self.num)

    def calc_sk_svk_smr1_smr2_sxp(self):
        '''
        sk: 中心部的水平差
        svk: 面的突出谷部深度
        smr1: 分离突出峰部与中心部的负载面积率
        smr2: 分离突出谷部与中心部的负载面积率
        sxp: 极点高度
        
        returns
        -------
        tuple-->sk,svk,smr1,smr2,sxp
        
        '''

        z_min, z_max = np.min(self.z), np.max(self.z)
        
        # 求给定高度c对应的smr值与smr=p的差值，用于二分法查找的被求解函数
        def fun_smr(c, p):
            return np.divide(self.z[self.z >= c].size, self.num) - p

        # 二分法查找smr=p对应的高度值c
        def binary_search(fun, p, start=z_min, end=z_max, precision=0.001):
            fa = fun(start, p)
            fb = fun(end, p)
            nums = 0
            if abs(fa) < 0.001:     # 判断区间边界值是否符合条件，近似零视为等于0
                return start
            if abs(fb) < 0.001:
                return end
            while abs(end - start) > precision and nums < 20:
                mid = (start + end) / 2.0
                fmid = fun(mid, p)
                if fmid == 0:
                    return mid
                if fa * fmid < 0:
                    end = mid
                if fmid * fb < 0:
                    start = mid
                nums += 1
            return start
        
        diff_c1_c2 = []     # 记录smr间距40%的高度差值
        for i in np.arange(0, 0.6, 0.01):
            c_1 = binary_search(fun_smr, i)
            c_2 = binary_search(fun_smr, i + 0.4)
            diff_c1_c2.append(c_1 - c_2)

        min_index = np.argmin(np.array(diff_c1_c2)) * 0.01  # 中心部分（最平缓）对应的端点smr值
        c1 = binary_search(fun_smr, min_index)
        c2 = binary_search(fun_smr, min_index + 0.4)
        # 直线方程 y=(y1-y2)/(x1-x2)*(x-x1)+y1 求sk
        c_smr1 = (c1 - c2) / (-0.4) * (0 - min_index) + c1
        c_smr2 = (c1 - c2) / (-0.4) * (1 - min_index) + c1
        
        sk = c_smr1 - c_smr2
        smr1 = self.calc_smr(c_smr1)[0]
        smr2 = self.calc_smr(c_smr2)[0]
        spk = z_max - c_smr1
        svk = c_smr2 - z_min
        #计算sxp，默认smr=2.5%与smr=50%对应的高度差
        sxp = binary_search(fun_smr, 0.025) - binary_search(fun_smr, 0.5)
        return sk, svk, smr1, smr2, sxp
                
    def calc_3d_args(self):
        '''计算面粗糙度参数'''

        sa = self.calc_sa()
        sv = self.calc_sv()
        sp = self.calc_sp()
        sz = self.calc_sz()
        sq = self.calc_sq()
        ssk = self.calc_ssk()
        sku = self.calc_sku()
        vvv, vvc, vmp, vmc = self.calc_vvv_vvc_vmp_vmc()
        sk, svk, smr1, smr2, sxp = self.calc_sk_svk_smr1_smr2_sxp()
        smr = self.calc_smr(0)[0]
        sdq = self.calc_sdq()
        sdr = self.calc_sdr()

        arg_3d=(
                    '%.3f'%sa,
                    '%.3f'%sv,
                    '%.3f'%sp,
                    '%.3f'%sz,
                    '%.3f'%sq,
                    '%.3f'%ssk,
                    '%.3f'%sku,
                    '%.3f'%vvv,
                    '%.3f'%vvc,
                    '%.3f'%vmp,
                    '%.3f'%vmc,
                    '%.3f'%sxp,
                    '%.3f'%sk,
                    '%.3f'%svk,
                    '%.3f'%smr,
                    '%.3f'%smr1,
                    '%.3f'%smr2,
                    '%.3f'%sdq,
                    '%.3f'%sdr)
        return arg_3d



if __name__ == '__main__':
    
    z_data = np.load('./temp/Z.npy')
    x_coord = 10	#任取一行
    
    surface_args_2d = Extract_2d(z_data, x_coord)
    print(surface_args_2d.calc_2d_args())

    surface_3d_args=Extract_3d(z_data)
    print(surface_3d_args.calc_3d_args())
